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Abstract 

We determine the stationary two-point correlation function of the one- 
dimensional KPZ equation through the scaling limit of a solvable micro- 
scopic model, the polynuclear growth model. The equivalence to a di- 
rected polymer problem with specific boundary conditions allows one to 
express the corresponding scaling function in terms of the solution to a 
Riemann-Hilbert problem related to the Painleve II equation. We solve 
these equations numerically with very high precision and compare our, up 
to numerical rounding exact, result with the prediction of Colaiori and 
Moore ^ obtained from the mode coupling approximation. 

1 Introduction 

In their well-known work [2] Kardar, Parisi, and Zhang argue that surface growth 
through random ballistic deposition can be modeled by a stochastic continuum 
equation, which in the case of a one- dimensional substrate reads 

dth = lX{d,hf + Ud^h + 7]. (1.1) 

Here h{x,t) is the height at time t at location x relative to a suitable reference 
line. ri{x,t) is space-time white noise of strength D, {'r]{x,t)ri{x' ,t')) = D6{x — 
x')6{t — t'), and models the randomness in deposition. z/9^/i is a not further 
detailed smoothening mechanism. The important insight of [2] is to observe that 
the growth velocity is nonlinear, in general, and is relevant for the large scale 
properties of the solution to (jl.lj) . To simplify, the growth velocity is expanded 
in the slope. The first two terms can be absorbed through a suitable choice 
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of coordinate frame. The quadratic nonlinearity in (ll.lj) is relevant and higher 
orders can be ignored, unless A = 0. 

The one-dimensional KPZ equation (jl.lll is regarded as exactly solved in 
the usual terminology. In fact, what can be obtained is the dynamic scaling 
exponent z = 3/2 jSlElllI- No other universal quantity has been computed 
exactly so far. In our contribution we will improve the situation and explain 
how to extract the scaling function for the stationary two-point function. A few 
other universal quantities can be computed as well. But they have been discussed 
already elsewhere 13 El E! • 

In 13] a mode-coupling equation for the two-point function is written down, 
in essence following the scheme from critical dynamics and kinetic theory. At 
the time only 2; = 3/2 and a few qualitative properties could be extracted from 
the mode-coupling equation. In jH] this equation is solved numerically. Such 
computations are repeated in Q with greatly improved precision and using a 
more convenient set of coordinates. Thus for the ID KPZ equation we are in 
the unique position of an exact solution and an accurate numerical solution to 
the mode-coupling equation with no adjustable parameters. As will be explained 
below, given the uncontrolled approximation, mode-coupling does surprisingly 
well. 

To attack (jl.lj) directly does not seem to be feasible, a situation which is 
rather similar to the one for two-dimensional models in equilibrium statistical 
mechanics. For example, the Ginzburg-Landau 0^-theory is given through the 
(formal) functional measure 



for the scalar field (p. (jl.2|) is not the proper starting point for computing the 
exact two-point scaling function at the critical coupling Qc- Rather one discretizes 
through the lattice and replaces the 0-field by Ising spins ±1. Then, following 
e.g. ini, the scaling function at and close to criticality can be obtained. By 
universality this scaling function is the one of (jl.2p . (While certainly true, to 
establish universality is difficult and carried out in a few cases only ^Ol-) In 
the same spirit we replace (jl.ip by a discrete model, where the most convenient 
choice seems to be the polynuclear growth (PNG) model. 

Before explaining the PNG model let us review the standard scaling the- 
ory for If the initial conditions h{x,0) of the KPZ equation are dis- 
tributed according to two-sided Brownian motion, then formally the distribution 
of h{x, t) — h{0, t) is again two-sided Brownian motion. Therefore it is natural to 
define the stationary time correlation 



where from the height difference the average displacement is subtracted. By 




(1.2) 




(1.3) 
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assumption 

C{x,0) = A\x\ (1.4) 

with roughness amphtude A = D/u to ensure stationarity in time, li z = 3/2, 
then C{x,t) scales as 

C{x,t) oct'^^^g {const ■ x/t'^^^), asx,t^cx) (1.5) 

with a universal scaling function g{y) having the asymptotics g{y) ^ Cq > for 
?/ — > and ^^d/) ~ Cooll/I for ||/| — * oo. In order to define g a.s a dimensionless 
function we fix the proportionality constants in p.Sj) as appropriate combinations 
of A and A, 

where the particular choice of numerical prefactors is in principle arbitrary. The 
factor 2~^/'^ in the denominator is chosen in order to conform with the convention 
for the GUE Tracy- Widom distribution The factor 2^/^ in the argument of 
the numerator differs from the convention used by Baik and Rains ^21 by a factor 
2 but conforms with the definition of the closely related Airy process j7j and has 
the further advantage to absorb a lot of prefactors in the equations defining g{y). 
Note however, that the exponents for the parameters A[^], i^fY']' ^[^] 
are fixed uniquely by dimensional reasoning. 

We remark that the slope dxh{x, t) is space-time stationary in the usual sense. 
For fixed t, x ^ dxh{x,t) is white noise with strength A. Since {dxh) = 0, the 
standard 2-point function is 

{d,h{Q, Q)dMx, t)) = \dlC{x, t). (1.7) 

This relation and the asymptotic behavior oi g, g{y)/\y\ — * 2 as y ^ oo, motivates 
the definition of a second scaling function, 

f{y) = l9"iy), (1.8) 

which by definition has integral normalized to one and which will be shown to be 
positive in the next section. 

In the sequel we will analyze the distribution function for the height differences 
in the stationary PNG model. As shown in jT2j, they can be represented in terms 
of certain orthogonal polynomials, which lead to recursion relations connected to 
the Painleve II differential equation ^3 • The asymptotic analysis is carried 
out in Our own contribution is twofold: (i) We observe that the stationary 
PNG model maps to a last passage percolation with boundaries [H]. (ii) The 
expressions in ^21 are given in terms of certain differential equations and the 
extraction of the scaling function g requires a careful numerical integration. This 
is one central point of our article. We will provide then plots of the structure 
function and give a comparison with the mode-coupling theory. 
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2 The polynuclear growth model 



The polynuclear growth (PNG) model is a model for layer-by-layer growth through 
deposition from the ambient atmosphere. The surface is parameterized by a time 
dependent integer- valued height function h{x,t), t G M, above a one-dimensional 
substrate, x G M. Thus the height function consists of terraces bordered by steps 
of unit height. The up-steps move to the left and the down-steps to the right 
with speed 1. Steps disappear upon collision. In addition to this deterministic 
dynamical rule new islands of unit height are nucleated randomly with space-time 
density 2 on top of already existing terraces. The corresponding stochastic pro- 
cess h{x, t) is well defined even in infinite volume (cf. for the closely related 
Hammersley particle process). 

Of interest to us here is the stationary growth process, which means that the 
slope dxh{x,t) = p{x,t) is stationary in space-time. One can think of p{x,t) as 
the density of a particle/antiparticle process. The particles are located at the 
up-steps and thus move with velocity —1, the antiparticles are located at the 
down-steps and move with velocity 1. Upon collision particle/antiparticle pairs 
annihilate. In addition, with space-time density 2, a particle/antiparticle pair is 
created with the particle moving to the left, the antiparticle to the right. To make 
p stationary, one prescribes at t = up-steps Poisson distributed with density p+ 
and down-steps independently Poisson distributed with intensity p_ such that 

p+p_ = 1. (2.1) 

This measure for steps is stationary under the PNG dynamics. The mean slope 
is given by 

M = p+-p_ = (9,/i(x,t)), (2.2) 

which is the only remaining free parameter. For fixed t, x ^—>- h{x, t) — h{0, t) is a 
(two-sided) random walk with rate p± for a jump from n to n ± 1. It has average 
u and variance p+ -|- p_, which implies for the roughness amplitude 

A{u) = V4 + m2. (2.3) 

For the growth velocity one obtains 

v{u) = {dth) = p_^ + p^ = V4 + m2. (2.4) 

Given p(x, t) the height h{x, t) is determined only up to a constant which we 
fix as /i(0, 0) = 0. To emphasize that only height differences count, /i(0, 0) is 
sometimes kept in the formulas. 

The stationary process with slope u transforms to the stationary process with 
slope through the Lorentz transformation 

x' = {l-cY^/\x-ct), t' = (l-c2)-i/2(t-cx), (2.5) 
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with the speed of "hght" equal to 1 and the velocity parameter c = — m/v4 + ??_ 
Thus it suffices to restrict ourselves to m = which we do from now on. In 
particular p+ = 1 = p_. ( ■ ) and E refer to the stationary density field at slope 
u = 0. 

The central objects are the height-height correlation 

C{x, t) = {{h{x, t) - h{0, 0) - 2tf) (2.6) 
and the closely related two-point function for the density, 

S{x,t) = {p{x,t)p{0,0)). (2.7) 

|92E(^(/i(x,t) -/i(0,0) -2t)') 

d^E(^p{x,t){h{x,t) - h{0,0) - 2t)) 

d^E (^p(0, t) (/i(0, t) - /i(-x, 0) - 2t) ) 
E(p(0,t)p(-x,0)) 

S{x,t). (2.8) 

The height correlation is convex, equivalently 

S{x,t)>0. (2.9) 

To prove this property we show that the structure function S{x, t) can be regarded 
as the transition probability for a second class particle starting at the origin. 
Its initial velocity is ±1 with probability |, as for the "ffist-class" up/down- 
steps. In contrast to an ordinary step the second class particle is never destroyed 
upon colliding with another step. Rather it eats up the step encountered and, 
by reversing its own direction of motion, continues along the trajectory of the 
absorbed step, cf. Figure Q Let p(x, t) be a given realization of the PNG process. 
The second class particle is added as 

p^^\x,0) = p{x,0) +aS{x), a = ±l. (2.10) 

p*-'^^(x, 0) evolves to p^'^\x, t) with nucleation events identical to the one for p{x, t). 
By construction, if Xt denotes the position of the second class particle at time t, 

p^''\x,t)- p{x,t) = a6ix-Xt). (2.11) 

Noting that by the Poisson property p'^'^^(x, 0) is given by p(x, 0) conditioned on 
the presence of either an up-step (a = +1) or down-step (cr = —1) at the origin. 



They are related as 

'^dlC{x,t) = 
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Figure 1: The trajectory of a second-class particle. 



we obtain 

< pt{x) 



(T=±l a=±l 

lim cr^ (p(a;, t) f piv, 0) dy = 

(T=±l -^-'5 



1 V- ,¥{jlsPix + y,t)dy = a'j'^^piy,0)dy = a} 

lim - y aa j 

'^°2^,.^^^ 26F{f_,piy,0)dy = a} 

iE(p(x, t)p(0, 0)) = t). (2.12) 

For arbitrary slope u the normalization of S{x,t) would be given by v{u) 



\/4 + and the mean oipt{x) evolves along the characteristics of the macroscopic 
evolution equation dtu = —dxv{u). Thus 

J S{x, t)dx = V4: + u"^, and J x S{x,t)dx = —tu. (2-13) 



3 The distribution functions for the height dif- 
ferences 

For the PNG model the distribution function for the height difference h{x, t) — 
h{0, 0) satisfies certain recursion relations, which are the tool for analyzing the 
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scaling limit when t oo and x = yt"^^^ with y = 0{1). The second moments 
yield C{x,t) and therefore by (j2.8p also S{x,t). 

Since the nucleation events are Poisson, h{x, t) — h{0,0) depends only on the 
events in the backward light cone {{x',t') G M^, < t' < t, \x — x'\ < t — t'} 
and the initial conditions at t = 0. Along the line {x' = t'} the down-steps are 
Poisson distributed with line density and correspondingly for the up-steps 
along the line {x' = —t'}. This property can be deduced from the uniqueness of 
the stationary state at given slope and the Lorentz invariance fl2.5|) in the limit 
c —>■ ±1. Thus h{x,t) — h{0,0) is determined by the nucleation events in the 
rectangle R^^t = {{x',t') G M?, \x'\ < \t'\, \x — x'\ < t — t'} together with the 
said boundary conditions. h{x,t) — /i(0,0) can be reexpressed as a directed last 
passage percolation according to the following rules: Inside R^^t there are Poisson 
points with density 2. Along the two lower edges of R^^t there are independently 
Poisson points with line density -\/2. A directed passage from (0, 0) to {x, t) 
is given through a directed path (polymer). It is a piecewise linear path in the 
plane, starts at (0, 0) and ends at {x, t), alters its direction only at Poisson points, 
and is time-like in the sense that for any two points {x',t'), {x",t") on the path 
one has \x' — x"\ < \t' — t"\. Note that, once the directed path leaves one of the 
lower edges to move into the bulk, it can never return. By definition the length 
of a directed path equals the number of Poisson points traversed. With these 
conventions 

h{x, t) — h{0, 0) = maximal length of a directed path from (0, 0) to {x, t). 

(3.14) 

We remark that in general there are several maximizing paths, their number 
presumably growing exponentially with t. 

Under the Lorentz transformation (j2.5|) the distribution for the height differ- 
ences (j3.14j) does not change. Therefore, we might as well transform R{x,t) to a 
square. By an additional overall scaling by V2 one arrives at av x v square, v = 
— x^, with bulk density 1 and the line densities a_ = a = {t — x){t + x) 
for the lower left, resp. = 1/a for the lower right edge. In this way we have 
recovered precisely the setting in IT^ with t replaced by v. Baik and Rains 
derive an explicit expression for the height distribution in terms of Toeplitz deter- 
minants, which can be further simplified by means of corresponding orthogonal 
polynomials. 

Let us state the result for the distribution function of h{x, t) — h{0, 0), 

F^^tin) = ¥{h{x,t)-h{0,0) <n} 

= Gnia)F{n)-Gn^iia)F{n-l). (3.15) 

For fixed v the functions g and F are given in terms of the monic polynomials 
TTniz) = z"' + 0{z^~^), which are pairwise orthogonal on the unit circle \z\ = 1 
with respect to the weig ht e^(^+^ '\ Their norm Nn is given by 

(7rn,Vrm) = Sn,mNn (3.16) 
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with (p, g) = § p{z)q{z i)e"*-^"'"^ ^\2T:iz) ^dz. One has 

n-l 



F(n) = e-^' JJiVfe, (3.17) 

fc=0 

where F(n) itself is the distribution function of the maximal length of a directed 
path in the case a+ = = a_. Thus lim^^oo F{n) = 1 and 

n 

G'„(a) = e-"(°+"")iV„ J]iV,7Vfc(-a)7rfc(-a-i) 

fc=0 

= e-''("+""') ((1 - n)7r„(-a)7r„,(-a-i) 

— Q;7r^(— a)7r„(— — a~"'^7r„(— a)7r^(— a^^) ) . (3.18) 



Defining the dual polynomials 7r*(2;) = z^T^n^z ^), the second equality in ()3.18j] 
is an easy consequence of the Christoffel-Darboux formula |16j . 

^ 7rfc(a)7rfc(fe) _ <(aX(&) - 7r„.(a)7r„(fe) 

fc=0 

valid for a, 5 G C, a6 7^ 1 and extended by I'Hospital's rule to ah = 1, and the 
trivial relation 

7i:{z)z-\:'{z-') + zn'^{z)n^{z~')=nn^{z)nn{z-') = n<(z)<(z-i). (3.20) 

Taking only the leading order of a in ()3.19j) one obtains the well-known rela- 
tions 

nn+l{z) = Z7ln{z) + Pn+inl{z), 

K+li^) = ZPn+lTTniz) +71*^{Z), (3.21) 

= iV„(l-p^+i) 

which are closed given Pn = vr„(0) for n > 0. For the particular weight function 
Q^iz+z ) Qj^g gg^j^ derive a nonlinear recursion relation for the PnS, 

Pn = (j)n+l+Pn-l){l-pl), (3.22) 

n 

with initial values po = 1, pi = h{2v) = {2ti)~^ ^^"^ e'^^e^^'^^^'^de is the 

modified Bessel function of order k and thus A^o = -^o(2'y). Eq. (j3.22|) is the 
discrete Painleve II equation. It has been derived in the context of orthogonal 
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polynomials for the first time in and later on more or less independently in 
[T71 UHl ESI Cni • The differential equations for 7r„,, vr*, 

7r„'(z) = {n/z + V/Z^ - Pn+lPnV/z)'Kn{z) + [pn+lV/z -pnV/z^)'nl{z) 
^n'(^) = {-Pn+lV/z+PnV)'Kn{z) + {-V +Pn+lPnV/z)'Kl^{z), (3.23) 

can be shown to hold by a tedious but straightforward induction, using ()3.2H) and 
()3.22|1 . They are implicitly derived in ^Hj, from which we learned their actual 
form. In [213] an integral expression is obtained for the derivative of orthogo- 
nal polynomials on the circle with respect to (up to some technical conditions) 
an arbitrary weight function. Specializing to the weig ht e^(^+^ resuhs m a 
differential-difference equation equivalent to ()3.23|1 . 

Of course, the mean of the probability distribution F^^tin) — F^^tin — 1) is 2t 
and its variance, the correlation function ()2.fi|l . is given by 

C{x,t) = J2 (2(^ - 2t) - l)F,,t(^). (3.24) 

n>0 

Thus to establish ()1.5p . one has to understand the scaling properties of the dis- 
tribution function F^^tin). Let us introduce the new variables s, y defined by 

n = 2v + v^/h, (3.25) 
X = v^^^y, (3.26) 

where v = \Jt'^ — x'^ is regarded as fixed when varying n and x. In ^2j the 
different scaling variable w = \y used, which leads to a string of factors of 2, 
avoided by our convention. Setting 

= -(-l)>n, (3.27) 

we rewrite fj3.22j) as 

(f - 2)Rn + 2Rl 

Rn+l — 2-R„ + Rn-1 = — • (3.28) 

i-Ri 

Under the scaling ()3.26|) . i?„ = v~^/'^u{y~^^^{n — 2v)) + C(t>~i), it becomes the 
Painleve II equation 

u'\s) =2u{sf + su{s), (3.29) 

in the limit v ^ oo. The starting value Rq = —1 is consistent with the left 
asymptotics of u{s) only if 

u{s) ~ — a/— s/2 as s — > — oo, (3.30) 

which singles out the Hastings-McLeod solution to ()3.29|) pi] . This particular 
solution will be denoted by u{s) and we conclude that 

u{s) = lim v^^^Ri2y^^i/3s], (3-31) 
9 



provided the limit exists (A complete proof is the main content of |221)- u{s) < 
and u has the right asymptotics 

u{s) ~ — Ai(s) as s ^ oo. (3.32) 

At this point we can derive the scaling limit for F[n). It is the GUE Tracy- Widom 
distribution function [TT] 

oo 

2 



^buE(s) = e-^(^), V{s) = - / v{x)dx, v{s) = {u{sf + s)u{sf - u'{ 



s 



(3.33) 

which appears already as the limiting height distribution for nonstationary curved 
selfsimilar growth [21]. Since v'(s) = u(s)^, one has v^^^ \ogN^2v+v^/^s] ~^ ^i^) 
F{[2v + v^/^s]) Fgue(s) as w ^ oo. 

Next we turn to the scaling limit for the orthogonal polynomials. For it to be 
nontrivial we set 

Pn{a) = e-''"<(-«), 

Qnia) = -e-^"(-l)"7r„(-«). (3.34) 
()3.26|1 implies a = 1 — ir^/^y + 0{v~'^^^). Setting n as in ()3.25j) . we claim that 
a{s,y) = Urn Pn{a), b{s,y) = Urn Qn{a). (3.35) 

If SO, the limit functions a, b satisfy the differential equations 

dstt = ub, 

dgb = ua — yb, (3.36) 

as a consequence of ()3.2H) . and 

dya = u^a — {u' + yu)b, 

dyb = (u — yu)a + {y'^ — s — u'^)b, (3.37) 

as a consequence of ()3.23|) . From ()3.2H) one immediately obtains 7r*(— 1) = 
(-l)"7r„(-l) = nLi(l - Rk)- One has the limit UT=ih - Rk) = e" since 
g-D +''A''o nfc=i ~ Rk)~^ has an interpretation as a probability distribution 
function [21]. Therefore the initial conditions to ()3.37|) are 

a(s,0) = -6(s,0) = e-^^'\ U{s) = - u{x)dx. (3.38) 

J s 

The scaling limit of Gn{<y) as defined in Eq. (j3.18|) is the function 



G{s,y) = / a{s',y)a{s',-y)ds' 

J — oo 

= «(-s, -y)dya{s, y) - b{s, -y)dyb{s, y), (3.39) 
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where the second equahty can be verified by differentiation with respect to s and 
using the identity 

a{s,y) = -bis,-y)e-y-'y, (3.40) 

itself being a direct consequence of ()H.H7|1 and (j3.38j) . Putting these pieces to- 
gether we obtain as scahng hmit for the distribution functions F^^tin), 

Fy{s) = ^ {G{s + y', y)FGVE{s + y^)) . (3.41) 

The shift in (j3.4H) by y"^ comes from the fact that Fy{s) is evaluated for constant 
t = V + \v^'^y + 0{v^^/^). 

In conclusion we arrive at the scaling function g{y) as defined in the Intro- 
duction. From (fO)|) . with A = ^ and A = 2 for the PNG model, and (PT^ we 
obtain 

g{y) = j s^dFyis). (3.42) 

As already mentioned, except for ()3.23j) . all our relations are derived in ^21 
Ej, and the existence of limits is proven with Riemann-Hilbert techniques. For 
completeness let us collect some more properties of a shown in 



a{s,y) - 
a{s,y) - 
a{{2yY/^x + y^y) - 

a{{-2yy/^x + y',y) - 
Therefore Fy{s) is asymptotically Gaussian and we recover g{y) ^ 2\y\ for large 

y- 



1, as s 


+00, 


0, as s 


—oo, 


1, as y 


+00, 






(27r)i/2 j 


— oo 



4 Numerical determination of the scaling func- 
tion 

The key object in determining the scaling functions g{y), f{y) is the Hastings- 
McLeod solution [21] to Painleve II, u{s), which is the unique solution to 

u" = 2u^ + su (4.1) 

with asymptotic boundary conditions (j3.3()j) and (j3.32|l . Tracy and Widom 
[m 123] integrate ()4.H) numerically with conventional differential equation solvers 
using the known asymptotics at s = ±oo. The precision achieved with this 
technique does not suffice for our purposes, since we need u{s) as starting values 
(|3.38p for the differential equations (|3.37p . We develop here a different method to 
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obtain u{s), in principle with arbitrary precision. Next the functions a{s,y) and 
b{s, y) have to be determined, which directly leads to values for the distribution 
functions Fy{s). They have to be further integrated with respect to s in order 
to obtain their variance, which is the desired scaling function g{y). The Taylor 
expansion method to be explained intrinsically produces not only function values 
at a point but also higher derivatives. Therefore we obtain f{y) not by numer- 
ically differentiating g{y) but rather by direct calculation via the knowledge of 



In a first step, to obtain reliable approximations to the Hastings-McLeod so- 
lution, we need to guess its initial data at some finite sq by using asymptotic 
expansions around ±00. Any initial data u{so) = Uq, u'{sq) = Ui give rise to a 
maximal solution, u{s) of (jOJ), which admits analytic continuation to a mero- 
morphic function on C. The only essential singularity for Painleve II solutions lies 
at 00. If u{s) is close to u{s) we can estimate the difference A(s) = u{s) —u{s) by 
linearizing (j4.ip around u{s). One obtains that A(s)/A(0) is of order exp{D{s)) 
if (A(s), A'(s)) is in the unstable subspace and of order exp{—D{s)) for the sta- 
ble subspace, where D{s) ^ — 1(— 2s)^/^ for s <^ and ^ f^^^^^ for ■§ 3> 0. Note 
that on the exponential scale we are looking at, derivation with respect to s 
leaves invariant the order. Therefore the exponential orders of A(s) and A'(s) 
are the same. Since generically initial values at sq have a component in both 
subspaces one obtains that A(s)/A(so) is of order exp(|D(s) — D{sq)\) in the 
range of validity of the linear approximation, |A(s)| <^ |w(s)|. 

It turns out that the left asymptotic expansion in (— optimally trun- 
cated at large negative sq, gives rise to initial values with A(so) only of order 
exp(— 1(— 2so)^/^). Thus control of the approximation always breaks down near 
s = 0. Approximations of the right asymptotics on the other hand allow a, in 
principle, arbitrary precision on any given finite interval. 

For s —>■ 00 the deviations of u{s) from the Airy function can be expanded in 
an alternating asymptotic power series with exponentially small prefactor. 





Ai(.) 



32vr3/2s7/4 




(4.2) 



The coefficients are oq = 1, ai = ^2 
recursion relation 



1493 
1152' 



. , and can be obtained via the 
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)(n 




for n > 



(4.3) 



with initial conditions a_i = a_2 = 0. 




(4.4) 



0<k<l<n 



are the coefficients in the asymptotic expansion of Ai(a;)^ and 



(6n- l)(6n-5) ^. 

Ai„ = ^ ^ ^Ai„„i, Aio = 1 

72n 



(4.5) 
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are the coefficients of the asymptotic expansion of the Airy function itself 



e 3* (—1) 



n 



^ n>0 ^3 ^ 

Empirically we observe that for sq ^ the optimal truncation in ()4.2|) is 
^ ^Sq^"^ leading to an exponentially improved (relative) precision 



Unght,n{^0) - U[So) 



UiSo 



exp(-|sf). (4.7) 



The linear perturbation argument for u{s) with initial values u{sq) = Mright,n(so), 
u'{so) = liright.nl'^o), ^ = [f-^o^^], IS now Valid for a slightly smaller interval than 
[— 2so, For example, by choosing the interval [— sq] and gluing u{s) 

at the boundaries smoothly to the optimally truncated asymptotic expansions, 
the maximal relative error of u{s) with respect to the Hastings-McLeod solution 
is of order e~^/^'^o . For our purpose it turns out that we do not need values 
with s < —20. On the other hand, to access large values of y in (|3.41|) . because 
of the shift in ()3.4H) . we need u{s) for large s with high precision. Uright,n('5) 
is numerically costly to evaluate for large s, so we finally choose sq = 100 and 
integrate in the interval [—20, 200]. This requires a maximal working precision of 
1500 digits and, given the integration of ()4.1|) is precise enough, u{s) and u{s) 
coincide in the ffist 1000 digits for s G [—20, 115] and still up to 50 digits at 
s = 200, where u{s) ~ 10"^"^^. In the sequel we drop the distinction between u{s) 
and its numerical approximation. The arithmetic computing is done partially 
with Mathematica® and for the computationally intensive tasks with the C++- 
based multiprecision package MPFUN++ |2H]- 

To solve initial value problems for ordinary differential equations highly so- 
phisticated iteration schemes are available, like Runge-Kutta, Adams-Bashford 
and multi-step methods. For arbitrary high (but fixed) precision results, all these 
methods become ineffective, since the step size is a decreasing function of the re- 
quired precision goal for the solution and tends to become ineffectively small. The 
only remaining choice is to Taylor expand the solution at a given point. The step 
size is limited by the radius of convergence only and the precision is controlled 
by the error made in truncating the Taylor series at some order fT7\ . 

u{s) is expanded at sq as 



u[s) = 

n>0 



Y^Ms-soT- (4.8) 



For the Painleve II equation the expansion coefficients u„ at sq are determined 

by uq = u{so), ui = u'{sq) and 

2Mn^ + SoM„ + M„-l 
Un+2 = 7 — Tw — T , (4.9) 

(n + 2)(n + 1) 
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where Un^ = Y^^=o'^n-juf are the expansion coefficients of u{s)^ at sq, Un^ = 
Un- We include the factorial into the expansion coefficients instead of taking the 
bare Taylor coefficients, in order to reduce the workload from multiplications by 
binomials when multiplying two expansions numerically. 

Numerically we find that the Hastings-McLeod solution does not have any 
pole in a strip |Im(s)| < 2.9. To have a safety margin we choose a step size one 
for the extrapolation of the expansion (j4.8p . 

We take the starting values m(so), u'{sq) from ()4.2|) at sq = 100. The coeffi- 
cients of the functions f/(s), V{s), see fj3.33p and ()3.38|) . when expanded around 
So, as in ()4.8|) . are given by 

(2) 

n + \ [n + l)[n + 1) 

and Vi = Uq—uI+SquI, leaving unspecified the yet unknown integration constants 
Uq = U{so) and Vq = V{so). The higher expansion coefficients of u, u', U, 
and V are independent of the values for Uq and Vq and are calculated with the 
recursion relations ()4.9|1 and ()4.10|1 . The values of u, u', U, and at s = sq ± 1 
are extrapolated, the expansion coefficients at sq ± 1 iterated. Then values are 
calculated at s = sq ± 2 by extrapolation, and so on. A posteriori we assign to 
U{sq) and ^(so) values, such that ?7(200) = = ^(200). The numerical errors 
from iterating ()4.9|) and from truncating ()4.8|) can be neglected compared to the 
uncertainty originating from the initial conditions. At the end of this ffist step 
we have at our disposal the expansion coefficients for u, u', U, V at the integers 
in the interval [—20, 200]. For the convenience of the interested reader let us just 
state the values at s = up to 50 digits, 

n(0) = -0.367061551548078427747792113175610961512192053613139, 

n'(0) = 0.295372105447550054557007047310237988227233798735629, 

U{0) = 0.336960697930551393597884426960964843885993886628226, 

V{0) = 0.0311059853063123536659591008775670005642241689547838, 

which might be used as starting values for a quick conventional integration of 
Painleve II to reproduce parts of our results with much less effort but also less 
precision. Tables can be found at [22] • 

The next step is to determine a{s,y), b{s,y) at sq G {— 20, . . . , 200} in the 
interval y G [—9,9] employing ()3.37j) and ()3.38j) . Setting 

a'{s,y)= ^ a^,n(s - So)"'(y - z/o)", 

r?i,n>0 

bis, y)= ^-."(^ - 'or(y - ^o)"' (4.11) 

m,n>0 
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()3.37|) becomes a recursion relation for the expansion coefficients, 

^ m 

0'm,n+l = ~~r ^0'm-k,n — {k + l)Uk+lbjn~k,n ~ '^fc^m-fc.n-l) ; 

n + 1 ^-^ ^ ' 
k=0 

^m,n+l [ T ( ^m,n— 2 ^m— l,n 

n + 1\ 

m 

m—k,n ''^fe'^m— fc,n— l) ) ; 



k=0 



(4.12) 



n > 0, allowing one to determine ao,n, ^o,n, n > upon the knowledge of ao,o, 
6o,o- We integrate along ±y with an extrapolation step size of |. From one 
obtains the recursions 

m + 1 ^-^ 

k=0 

^ m 

&m+l,n = -( — bm,n-l + y Ukam-k,nl ■ (4-13) 

m + 1 V ^-^ / 



fc=0 



The expansion coefficients Gm,n of G{s,y) at (scyo); are determined from (j3.39|) 



as 



Gm,,n = (n + 1) (a„ „am,,n+l " 6m^„fcm,n+l) (4-14) 

where a~„, are the corresponding expansion coefficients of a and b at 

(■So, -Z/o)- 

To finally determine g{y) and its derivatives we write 



soez m>l 

Cm,„ are the expansion coefficients of {s,y) ^ J^^{r - y'^f£2{G{r,y)FGu^{r))dr 
at (so,?/o), 

Cm,n = (m - 2)(GF)„_i,„ - 2m(GF)„,„„2 + ^^5i^(GF)„+i,„_4. (4.16) 

Here 

(4.17) 

fc=0 

are the expansion coefficients of G{s,y)FQ\j-p^{s) and F„ = — Yl'k=i ^^kFn-k are 
the expansion coefficients of -Fgue- Numerically the sum over sq in ()4.15p is 
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truncated to values inside [—15,200], since outside contributions turn out to be 
negligible at the chosen precision goal. After accomplishing this program we keep 
values for g{y) aty e [-9, 9] and for ^("^ (?/) , n = 0, . . . , 4, at e |Zn [-9, 9] 

with an accuracy of about 100 digits (a table in ASCII format is available online 
at [2ni)- For interpolating these values we used the Interpolation-function of 
the Mathematica® package yielding best results due to the high precision data 
with an interpolation order of 57. 

5 Discussion of the scaling function 

There have been numerous attempts to approximately determine g{-) |^ 1301 1^ 
IH1E2]- For historical reasons a different scaling function, -F(-), is analyzed in some 
of these works. The relation to our scaling function g{-) is 

F{0 = i^/2f/'g{i2e)-'/'), resp. giy) = 2y F{l/i2'/Y^')). (5.18) 

Note that by ()1.6|) the large y behavior of g is fixed by definition as g{y) ~ 
2\y\. The special value ^(0) = 1.1503944782594709729961 is the Baik-Rains 
constant f2llE]- In the literature the universal amplitude ratio Rg = 2~'^^^g(0) = 
0.7247031092 and the universal coupling constant g* = ^(O)-^/^ ^ 0.810456700 
have been investigated. Approximate values have been determined by means 
of Monte-Carlo simulations for the single step model jS^, numerically within 
a mode- coupling approximation [HIH |Hl P, and even experimentally for slowly 
combusting paper [^31 yielding estimates for g{0) within reasonable ranges around 
the (numerically) exact value indicated. 

In Figure El the scaling function f{y) = \g"{y) is shown as determined by the 
multiprecision expansion method explained in the previous section. We estimate 
its large y asymptotics as 

log/(2/) ^ -c\y\^ + o{\y\) for y oo. (5.19) 

The cubic behavior is very robust and numerical fits yield about 2.996-2.998 quite 
independently of the assumed nature of the finite size corrections. The prefactor 
c = 0.295(5) has a relatively high uncertainty because of the unknown subleading 
corrections. Even though unaccessible in nature we estimate the error term, as 
indicated in ()5.19|) . to be sublinear or even only logarithmic from the numerical 
data. Possibly, the exact asymptotic behavior could be extracted from a refined 
asymptotic analysis of the Riemann-Hilbert problem. 

Colaiori and Moore ^ tackled the same scaling function by completely 
different means. Starting from the continuum version of the KPZ equation they 
numerically solved the corresponding mode-coupling equation |31 |H] , which con- 
tains an uncontrolled approximation, since diagrams which would renormalize the 
three-point vertex coupling are neglected. Nevertheless a qualitative comparison 
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Figure 2: The scaling function f{y) versus y in a semilog arithmic plot. The dotted 
line exp(— 0.295|y|^) is drawn as a guide to the eye for the large y asymptotics of 
/• 

of their result with the exact scahng function f{y) shows reasonable similarity, 
cf. Figure El Both functions are normalized to integral 1 by definition. The mode 
coupling solution oscillates around for \y\ > 3, whereas f{y)>0 for the exact 
solution. We do not know whether this is a numerical artifact or an inherent prop- 
erty of the mode-coupling equation. On the other hand, the second moments are 
reasonably close together, 0.510523 for f{y), and 0.4638 for the mode-coupling 
approximation. So is the value of the Baik-Rains constant g{0) = 2 J \y\f{y)dy 
for which mode-coupling predicts the value 1.1137. 

From the solution to the mode-coupling equations one does not directly obtain 
f{y), but rather its Fourier transform. The function G{t) from P is defined 
through 

C?(p/2/27/2) ^ ^ 2 / cos{ky)fiy)dy. (5.20) 

Moore and Colaiori predict a stretched exponential decay of G{t) as oc exp(— c|rp/^) 
[S3] and numerically find a superimposed oscillatory behavior on the scale |rp/'^ 
p. In Figure m f{k) is plotted as obtained by a numerical Fourier transform of 
f{y)- Indeed it exhibits an oscillatory behavior as can be seen in Figure El where 
the modulus of f{k) is shown on a semilogarithmic scale. The dotted line in the 
plot is the modulus of the function 

10.9A;-9/^ sin(iA:3/2 _ i_gs7)e-^'''^\ (5.21) 

shifted by a factor of 1000 for visibility, which fits f{k) very well in phase and 
amplitude for k ^ 15. This behavior is not in accordance with the results of 
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Figure 3: The exact scaling function f{y) compared to the mode coupling result 
of Colaiori and Moore ^(dotted line). Both functions are even. 



Colaiori and Moore, since the oscillations and the exponential decay of G{t) for 
the exact solution are apparently on the scale r and not r^/^. 

Note that f{k) is the scaling function for the intermediate structure function 

S{k,t) = j dxe'^^'Six.t) ~ 2/(^2/3^;). (5.22) 

By Fourier transforming with respect to t we determine the dynamical structure 
function, 

S{k,uj) = j cixdte'(^"+"*)^(x,t) - 2A;-3/2 j (uj/k^^^), (5.23) 

where 

/ (r) = I ds e^"/(.2/3) = 2 rfl/ T''L'{y/T'/')f{y) (5.24) 

and L has the convenient representation 

L{k) = 2 • 32/3Ai(-3-^/3/€2) sin(2KV27). (5.25) 
The correlation function ()2.6|) in Fourier space is given by 

C{k, uj) = 2k-^S{k, uj) ~ C^P^(A;, uj) = Ak-''^ f (00/^/^), (5.26) 

describing the asymptotic behavior at A;, u; = 0. Note that C{k,uj) > by defi- 
nition, since {hk^ujhk' ,uj') = 5k~k'5ui~u)'C{k,uj) for {k,uj) 7^ (0,0). The anomalous 



18 




Figure 4: The Fourier transform f{k) of the scaling function f{y). 



scaling behavior in real space is reflected by the exponents for the divergence of 
C^^'^{k,uj) at k = LJ = 0. In the linear case, the Edwards- Wilkinson equation 
A = in (jl.lj) . one easily obtains 

C^'^{k,u)= (5.27) 

A 3d-plot of C^^'^ik, oj) is shown in Figure IHl Its striking features are the smooth 
behavior away from fc, a; = 0, especially on the lines where /c = and u; = and 
the two symmetric maxima of t— > C^^^(/c, cj) for constant uj. Our numerical 
data yield for the singular behavior at A; = 0, u; = 0, 

= cc.-^/=^(2.10565(l) + 0.85(1) + 

= A;"^/2(19.4443(1)- 52.5281(1) cu2A;"=^ + 0(cj^A;"^)). (5.28) 



6 Conclusions and Outlook 

For systems close to equilibrium many properties valid in generality rely on de- 
tailed balance, amongst them in particular the link between correlation and re- 
sponse functions. The KPZ equation does not satisfy detailed balance, since the 
growth is directed. However, it has been speculated that in 1 + 1 dimensions de- 
tailed balance is recovered in the scaling regime. With our exact scaling function 
at hand, such a claim can be tested. 

Detailed balance implies that the eigenvalues of the generator in the mas- 
ter equation lie on the negative real axis. Thus autocorrelations in the form 
(X(t)X(O)) can be written as the Laplace transform of a positive measure. The 
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Figure 5: The modulus of f{k) on a semilog arithmic scale. The dotted line is a 
heuristic fit, shifted by a factor 1000 for visibility. 

structure function S{k,t) at fixed k is such an autocorrelation. Using the scahng 
form (|5.22|) detailed balance would imply 



with v{dX) > 0. In particular S{k,t) > 0. From ()5.21|) we know that S{k,t) 
oscillates around zero. Definitely, at \k\ ~ 5 there is a negative dip, cf. Fig. |3] 
Thus ()6.H) cannot be correct. 

The Bethe ansatz [3211201 indicates that, for large system size, the density of 
states is concentrated on an arc touching 0. If so, the integration in (j6.ip would 
have to be replaced by a corresponding line integral. It is not clear to us how to 
extract from the numerical knowledge of S{k,t) such a representation. 

Our main result is the exact scaling function /, see Figure El for the two-point 
function of the stationary KPZ equation in 1 + 1 dimensions. "Exact" must be 
qualified in two respects. Firstly / is given indirectly through the solution of 
certain differential equations, which can be solved numerically only with con- 
siderable effort. The errors are well controlled, however. Secondly, we rely on 
universality, in the sense that the scaling function is derived for the PNG model, 
which is one rather particular model within the KPZ universality class. Of course, 
it would be most welcome to establish the scaling limit also for other models in 
this class. 

The KPZ equation ()1.1|) is a two-dimensional field theory and, in spirit, be- 
longs to the same family as two-dimensional models of equilibrium statistical 
mechanics, one-dimensional quantum spin chains, and other (1 + l)-dimensional 




(6.1) 
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Figure 6: The correlation function C^^^{k,uj) in Fourier space. 



quantum field theories at zero temperature. While in the latter cases, there 
are a number of models for which the two-point function can be computed, in 
the dynamical context such solutions are scarce. In addition, the KPZ equation 
does not satisfy the condition of detailed balance. Such nonreversible models are 
known to be difficult and we believe that the PNG model is the first one in the 
list of exact solutions, disregarding noninteracting field theories. 

For the nonstationary KPZ equation with a macroscopic profile of nonzero 
curvature the analogue of Fq is the Tracy- Widom distribution function. In that 
case the full statistics of x i— h{x,t) for large but fixed t is available [7]. It is 
conceivable that an extension of the techniques used there also admits a more 
detailed study of, say, the joint distribution of h{x, t) — h{0, 0), h{x', t) — h{0, 0). 
On the other hand the joint distribution of h{0, t) — h{0, 0), h{0, t') — h{0, 0) does 
not seem to be accessible. In the representation through the directed polymers 
it means that space-like points, even several of them, can be handled, whereas 
time-like points remain a challenge. 
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